Accepted to ApJ: Mar 27, 2012 

Preprint typeset using I^'T^]X style eniulatoapj v. 5/2/11 



THE EVIL-MC MODEL FOR ELLIPSOIDAL VARIATIONS OF PLANET-HOSTING STARS AND 

APPLICATIONS TO THE HAT-P-7 SYSTEM 



Brian K. Jackson^ 

Carnegie Institution for Science, Washington, DC 20015, USA 



(N 

o 

(N 

u 

< 
0^ 



Oh! 
O 



(N 
> 

O 

o 

rn 

O 



1 



NiKOLE K. Lewis 

Lunar and Planetary Laboratory, University of Arizona, Kuiper Space Sciences Building, Tucson AZ 85721, USA 

Jason W. Barnes 

Department of Physics, University of Idaho, Engineering-Physics Building, Moscow, ID 83844, USA 

L. Drake Deming 

Department of Astronomy, University of Maryland, College Park, Maryland 20742, USA 

Adam P. Showman 

Lunar and Planetary Laboratory, University of Arizona, Kuiper Space Sciences Building, Tucson AZ 85721, USA 

AND 

Jonathan J. Fortney 

Department of Astronomy and Astrophysics, UCO/Lick Observatory, University of California, Santa Cruz, CA 95064, USA 

Accepted to ApJ: Mar 27, 2012 

ABSTRACT 

We present a new model for Ellipsoidal Variations Induced by a Low-Mass Companion, the EVIL- 
MC modeJ3- We employ several approximations appropriate for planetary systems to substantially 
increase the computational efficiency of our model relative to more general ellipsoidal variation models 
and improve upon the accuracy of simpler models. This new approach gives us a unique ability to 
rapidly and accurately determine planetary system parameters. We use the EVIL-MC model to 
analyze Kepler Quarter 0-2 (QO-2) observations of the HAT-P-7 system, an F-type star orbited by 
a ~ Jupiter-mass companion. Our analysis corroborates previous estimates of the planet-star mass 
ratio q = (1.10 ± 0.06) x 10^'^, and we have revised the planet's dayside brightness temperature to 



2680 



+ 10 
'20 



K. We also find a large difference between the day- and nightside planetary flux, with little 
nightside emission. Preliminary dynamical-|-radiative modeling of the atmosphere indicates this result 
is qualitatively consistent with high altitude absorption of stellar heating. Similar analyses of Kepler 
and CoRoT photometry of other planets using EVIL-MC will play a key role in providing constraints 
on the properties of many extrasolar systems, especially given the limited resources for follow-up and 
characterization of these systems. However, as we highlight, there are important degeneracies between 
the contributions from ellipsoidal variations and planetary emission and reflection. Consequently, for 
many of the hottest and brightest Kepler and CoRoT planets, accurate estimates of the planetary 
emission and reflection, diagnostic of atmospheric heat budgets, will require accurate modeling of the 
photometric contribution from the stellar ellipsoidal variation. 

Subject headings: Planets and satellites: fundamental parameters - Planets and satellites: individual: 
HAT-P-7 



1. INTRODUCTION 

The Kepler and CoRoT missions have begun a new 
chapter in time-domain astronomy. Among other re- 
sults, the phenomenal photometric stabilities, long obser- 
vational baselines, and high duty cycles of these missions 
will provide a vast harvest of new exoplanets. Already, 
the Kepler mission has found 25 planets that have been 
confirmed and an additional 1,235 planetary candidates 
porucki et al.ll2011f ). The stability of the Kepler and 
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CoRoT photometry also allows access to astrophysical 
signals with amplitudes too small to have been detected 
previously. 

Too small to have been observed before, the photo- 
metric signal of tidal distortion of a star by a close-in 
planet can now be measured using Kepler and CoRoT 
data. This signal is usually referred to as an ellip- 
soidal variation since a tidally dis torted body takes 
an approximately ell ipsoidal shape. IDrakd ()20Q3D and 
ILoeb fc Gaudil ()2003f ) initially suggested that the Kepler 
mission might observe ellipsoidal variations for many of 
its targets, and the latter study estimated amplitudes as 
large as 100 parts per million (ppm) for very short-period 



hot Jupiters. 

The amphtude of the ellipsoidal variation depends on 
several key system parameters, including the ratio of the 
stellar radius to the orbital semi-major axis a, the planet- 
star mass ratio q, and the sine of the orbital inclination 
sini. If a planet transits its host star, the transit light 
curve allows accurate determination of a and sini. If el- 
lipsoidal variations can also be measured for the system, 
the mass ratio itself can be estimated. Hence, with an es- 
timate of the stellar mass, ellipsoidal variations can help 
confirm the planetar y nature of a transiting companion 
(|ShDorer et al.llMll) . 

Kepler and CoRoT observations can also provide con- 
straints on the planetary emission and reflection, which 
can elucidate a planet's atmospheric properties. How- 
ever, in visible wavelengths monitored by the missions, 
the contrast between light emitted from a hot close-in 
planet and a host star is much smaller than it is in 
the infrared. Consequently, determination of a close-in 
planet's emitted and refiected flux requires accounting 
for the stellar ellipsoidal variation. 

Analysis of ellipsoidal variations and eclipses has a 
long history for close binary stars, where it provides 
a wealth of information regarding stellar masses, lumi- 
nosities, and internal structures, among other properties 
(|Kopallll959D . The effects of tides in such systems (tidal 
distortions, thermal perturbations, etc.) can be dra- 
matic, but generations of binary star astronomers were 
hampered by limited, semi-analytic models. The com- 
putational power required for highly accurate numerical 
models was o nly developed in the last few decades (e.g. 
iWilsonI [19941) . However, for planetary systems, tidal ef- 
fects are much smaller, owing to the small planet-star 
mass ratio {q < 0.01), and so they give rise to much 
smaller (and possibly less complex) ellipsoidal variations. 
Consequently, a model for ellipsoidal variations in plan- 
etary systems can be greatly simplifled relative to more 
general models appropriate for binary stars. 

In this paper, we present a new model for ellip- 
soidal variations in planetary systems - the Ellipsoidal 
Variations Induced by a Low-Mass Companion (EVIL- 
MC) model. We incorporate many approximations ap- 
propriate for planetary systems, which allow our model 
to be computationally efficient. Our model uses an al- 
ternative appr oach to other rece ntly applied or devel- 
oped models. ^W elsh et al.l (|2010 [) discovered ellipsoidal 
variations in Kepler observations of the HAT-P-7 sys- 
tem and analyzed them us ing the binary star ELC code 
(jOrosz fc Hauschildti[2000l ). That code is state-of-the- 
art but requires considerable computational resources to 
model the v ery small planet-induced ellipsoidal variation. 
iMazeh fcF aiglcr (2010) proposed a semi-analytic model 
involving a Fourier expansion of photometric signals in- 
duced by the pr esence of a planet, in clud ing the ellip- 
soidal variation. iShporer et al.l (|201lD and IMazeh et al.l 
()2011|) applied that model to photometric variations ob- 
served for the KOI-13.01 system. The simplicity of this 
model allows rapid analysis of data, but the relationships 
between the Fourier coefficients and the system parame- 
ters are not all accurately determined, limiting the ability 
of this model to determine the param eters. Specifically, 
in their analysis, IShporer et al.l (1201 1[ ) did not report a 
planet-star mass ratio based on the ellipsoidal variation. 
Determining that relationship requires a model that ac- 



counts in detail for tidal effects. 

For this paper, we tailor our model to Kepler observa- 
tions of the HAT-P-7 system and constrain the planet's 
mass, phase function, whence we derive constraints on 
the atmospheric brightness temperatures. We also high- 
light the importance of considering the ellipsoidal varia- 
tions when modeling planetary emission and reflection in 
visible wavelengths. In Section[2l we describe our model, 
derive the relevant equations, and compare our model to 
others. In Section[3l we describe the Kepler observations 
of the HAT-P-7 system and how we conditioned the data 
for analysis. In Section JH we apply the EVIL-MC model 
to these data. Finally, in Section [SJ we discuss implica- 
tions of our results and future work. 

2. MODEL DESCRIPTION 

In this section, we flrst describe the approximations 
made in our model and derive the relevant equations. 
Then, we compare our model to others. 

2.1. Model Approximations 

To model tidal distortion of stars hosting close-in plan- 
ets, we make several approximations: 

1. We treat the planet and star as point masses to 
determine their gravitational flelds, ignoring the 
contribution of the asymmetric mass distributions 
arising from tida l distortions , which is negligible 
for our purposes (|Claretll2000[ ). 

2. We employ the equilibrium tide approximation 
for modeling the stellar shape, in which the stel- 
lar surface lies along a gravitational isopotential 
(jWilson fc Sofial 119761 ). We neglect tidal dissi- 
pation or more complex hydrodynamic motions 
within the star that m ay be observable in some 
cases (jPfahl et al.l[200l . 

3. We assume that the planet's orbit is circular. Al- 
though the orbits of many tran siting planets are ec - 
centric (notably HD 80606 b - iWinn et aTll2009b( ). 
the majority are circular or nearly so. 

4. We assume that the planet-star mass ratio q is 
small and that the host star rotates as a solid body 
with a centrifugal acceleration that is small com- 
pared to the surface gravity. As a consequence, the 
departure of the stellar shape from sphericity due 
to tides and rotation is assumed small (tens of ppm 
for the HAT-P-7 system or - 10 km). 

5. We neglect Dopplcr effects from the velocity of tidal 
motions within the star. The rotational and tidal 
motions of the stellar su rface may contribu te to the 
Doppler flux variations ([Arras et al.ll20l3) but are 
probably negligible for broadband Kepler observa- 
tions. 

We do NOT assume that the stellar rotation is syn- 
chronous with the orbit or that the stellar obliquity is 
zero. Alignment and synchronization of stellar rota- 
tion with the orbit is common among close binary stars, 
likely as a result of large tidal torques (|Kopallll959D , but 
not for transiting planets. For the majority of planet- 
hosting stars for which it can be determined, the stellar 



rotatio n period is much lon ger than the planet's orbital 
period () Jackson et al.ll2008[). Also, observat ions of the 
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Rossiter-McLaughlin effect (jWinn et al.ll201Cl[ ). detection 
of the crossiiiK of s tar spots by transiting planets (e.g. 
iDeming et al.ll201ll ). and contrib utions of gravity d ark- 
ening to transit light curves (e.g. iBarnes et al.ll201ll ) all 
show that many planets have orbits that are strongly 
inclined relative to their host stars' equators. 

We also do NOT make the usual assumption in tidal 
modeling that the orbital separation is much larger than 
the physical radius of the star, i.e. we do NOT assume 
the star is an ellipsoid. Typically, the potential of a body 
inducing tidal distortion (and consequently the shape of 
the tidally distorted body) is expanded in the ratio of the 
ti dally-distorted body's radi us to the orbital separation 
a ([Murray fc Dermottlll999l ). For planets very close to 
their host star - the planets for which ellipsoidal variation 
will be the most pronounced - higher-order terms may 
contribute to the di stortion non-neg ligibly. For the HAT- 
P-7 system, a '^ 4 (jPal et al.l 120081 ). so depa rture of the 
stellar shape from an ellipsoid is significant. iPfahl et al.l 
(|2008| ) showed that assuming the star is an ellipsoid may 
not be sufficiently accurate for very close-in exoplanets, 
resulting in erroneous estimates of the system parame- 
ters. We discuss this point more in Section [2.31 

Currently our model does NOT include occultation of 
either the star (during transit) or the planet (during 
eclipse) . In fitting our model to data, we mask out these 
phases. A future version of our model will include these 
phases, but an initial analysis suggests that, for typical 
planetary systems, the tidal distortion of the host star 
negligibly modifies the transi t light curve. Thus, s tan- 
dard light curve models (e.s.. iMandel fc Agoill2002[ ) are 
probably sufficiently accurate. 

2.2. Model Equations 

We use a coordinate system centered on the star, as 
illustrated in Figure [TJ (In our notation, a vector Q 

has magnitude Q and is parallel to Q, a vector with 
unit length.) As measured in this frame (and subject to 
the approximations above), the gravitational potential 
on the stellar surface U is: 
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^i?* cos i; + -ujIrI (l - cos^ A) 
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where G is Newton's gravitational constant. Mi, the stel- 
lar mass, i?* the distance from the stellar center to its 
photosphere (which is not constant), Mp the planet's 
mass, A the orbital semi-major axis, w* the stellar ro- 
tation rate, cos?/' = R^, • A, and cos A = R* • w^,. The 
stellar rotation axis is w*, and the planet's position vec- 
tor is A. 

The first term in Equation[T]represents the star's gravi- 
tational potential, and the second term the planet's grav- 
itational potential. If the second term in Equation[l]were 
expanded as a Taylor series in Ri,/A, the first-order term 
(proportional to cos'i/') would correspond to the force 
constant throughout the star that keeps it in orbit about 
the system barycenter. Since this force is constant, it 
does not contribute to the tidal distortion, and so we in- 
clude the third term in Equation [T] to remove it. The 




Fig. 1. — Definition of model geometry. The large, light circle 
represents the star, and the small, dark circle represents the planet. 
The coordinate system (X, Y , Z) is centered on the star, Z points 
toward the observer, X points along the orbital line of nodes, and 
Y is in the plane of the sky and perpendicular to X. R,* points 
somewhere on the stellar surface, A points to the planet in its orbit 
(a portion of which is illustrated), and w* is the stellar rotation 
vector. Ro points somewhere on a stellar surface at a right angle 
to both uii, and A. The relevant angles are iji, the angle between 
Rj. and A, and A, the angle between R* and u)^. 



fourth term in Equation [T] represents a potential corre- 
sponding to the centrifugal acceleration due to the stellar 
rotation. 

At points on the stellar surface where R* _L A and 
R* _L w*, we take R* = Rq. To clarify this definition, 
consider the case of a planet crossing the exact center 
of the disk of a star with a rotation vector pointing ex- 
actly at the observer. At this instant, A || o)^ || Z, and 
i?o would be the usual stellar radius that goes into de- 
termining the transit depth. In the general case, the 
relationship between the transit depth and Rq is more 
complicated. However, the tidal distortion has a negli- 
gible effect on the transit light curve, so we can safely 
consider i?o as the usual radius that goes into determin- 
ing the transit depth. 



We normalize U by ( ^^f^ ) , giving $: 



$ = [/ 



\GM^J R ' [a? - 2aR cos ijj + R^ 



,1/2 



n 1 / ?^ /t*^ 

-\Rcosi: + - — ^ (1 + q)(l- cos^ A) (2) 

where R = i?*/i?o, 9 = Mp/Mi,, a = j4/i?o, and w = 
ijJi,ln, with n as the orbital mean motion. 

Normalized to Rq, the stellar radius i? = 1 + (5_R, a 
function of cos-i/j and cos A. Per our assumptions, the 
surface of the star corresponds to an isopotential contour, 
i.e. $ = const. We take the constant to be the potential 
$0 at Ro (where cosV' = cos A — 0): 



$0 = 1 



1 



(1/2 2 a- 



0-^(1 + 9). 



(3) 



The departure from sphericity (5R is assumed small, and 
we can expand $: 



$ 
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1 + bR (a2 _ 2a(l + 8R) cosiP + (1 + 6R)^)^^^ 

±{l + SR)cos^ + l- '^'^^\^^^\ l + q)(l-cos^X) 
a'' Z a-^ 



(1 - SR) + 
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2a3 
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1'^^ /1 2 \\ 
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(a2 + 1) 
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1/2 + 2 a3 ■ 
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In Equations S] and [SJ we have dropped 2nd-order 
terms. We set Equation 2] equal to [5] and solve for 5R: 



5R = q[ [a^ - 2a cos V + 1] 

^^ 2 ^ 

— COS A. 



-1/2 



1-1/2 



C0S1/' 



(6) 



Gravity darkening of the stellar surface also con- 
tributes to the photometric variation. Briefly, the 
planet's tidal gravity perturbs the balance of forces (pres- 
sure, the star's own gravity, radiation, etc.) within 
the stellar atmosphere and results in a small (few 0.1 
K) decrease in the effective temperature and bright- 
ness at points on the stellar surface nearest the planet 
(|yon Zeipel 1924!). Theoretical considerations motivate a 
parameterizatio n of the gravity darkening involving the 
surface gravity ([Wilson fc Deviniievlll971fi . The gravity 
vector on the stellar surface is given by 



g = — ^3^R 



GMp (A - R,) 



GAL 



^A 



■UJ_ 



Rl '"* ' {A^^2R^AcosiP + Rlf^^ A^ 
^i?, (r, -lD.cosA). (7) 

The last term in Equation [71 representing the centrifugal 
acceleration, is usually written as — w* x (uji, x R*). Using 



cross-product identities, this expression can be re- written 



R*(w* • w*) - ^^(w* • R^) = a;^i?*(R* - w^cosA). 
We normalize g by ( '^^f* ) , giving F: 



r = g 
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1 q ( aA — i?R^, 

rR* 



\GMJ i?2"* ' (a2_2ai?cos7/> + i?2)3/2 

R 
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( a A — R* 



(a2-2acosV^-f 1)^/^ 
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a^ a-^ 

= -R, + 6r (8) 

where ST represents all the gravitational accelerations 
other than the zeroth-order stellar gravity. The magni- 
tude of r is 



r = (r • r)^/^ 
= ([-R, + (5r]-[-R, + jr] 
w 1 - R* • 5r. 



1/2 



(9) 



Likewise, at Rq, the gravity vector is Fg « 1 — Ro-^Fq. 
Using these expressions, the effective temperature at T 
on the stellar surface is parameterized as 



r = rW — 1 ~T, (^l-|-/3[Ro-<5Fo-R. •(JFlj (10) 

where T is the temperature at R^, and /3 the grav- 
ity darkening exponent. T^, is the effective temperature 
at Rq. For our analysis, we are only interested in the 
fractional variation in the stellar brightness, and surface 
brightness variations are linear in the small difference 
in temperature between Rq and any other point on the 
surface. Moreover, tides raised by planets have a negli- 
gible effect on the determination of the stellar effective 
temperature from observation, and the usual distinctions 
between a sta r's polar and mean effective temperatures 
([Wilsonl [19791 ) are unimportant here. Consequently, we 
take T-t to be both the mean effective temperature and 
the temperature at Rq. 

To model the limb-darkening of the stellar disk, we 
calculate the projection of the normalized gravity vec- 
tor onto the line of sight, /z = F • Z. We use this to 
determine the li mb-darkening profile Iju) assuming a 
quadratic profile ([Mandel fc Ago]|l2002[ ): 



/(/i)//(l) = l-7i(l-M)~72(l-/i)' 



(11) 



where 7^ are the limb-darkening coefficients. The model 
DOES allow for other profiles, though. 

Our model also includes the photometric effects of the 
stellar refiex velocity vz, ref e rred t o as "Doppler fiux 
variations" in iLoeb fc Gaudil ([2003D . These variations 
come in at the first order in the ratio of vz to the speed 
of light and include several effects convolved together: 



(1) transformation of the energy-momentum four- vector 
from the frame co-moving w ith the star to the observer 's 
frame (Equation 4.93 from iRvbicki fc Lightmanlll979D . 

(2) reduction in the apparent size of the star as it re- 
cedes from the observer (Equation 4.95 from ibid.), (3) 
increased travel time for the steUar photons as the star 
recedes from the observer (see discussion point 2 above 
Equation 4.97 in ibid.), and (4) Doppler shifting of the 
steUar flux measured within the observational bandpass. 
Together, effects (l)-(3) increase the apparent stellar flux 
as the star approaches the observer. The peak in emis- 
sion for HAT-P-7 occurs blueward of the Kepler band- 
pass, so the accompanying blue-shift of the stellar flux 
(effect 4) reduces the apparent flux. However, taken al- 
together, the Doppler flux variations cause HAT-P-7 to 
brighten as it approaches and darken as it recedes. 

To first order in q, the star's line-of-sight velocity is 



vz 



- {q sin i)nA sin(27r0) 



27rGA'f* 
P 
Kz sin(27r 
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(qsini) sin(27r(5 



(12) 



where A is the orbital semi-major axis (NOT normal- 
ized to Ro), 4> is the orbital phase (= at mid-transit), 
P is the orbital period, and Kz is the amplitude of the 
projected reflex velocity of the star. (Not e that we have 
chose n the opposite sign convention from lLoeb fc GaudU 
l2003t positive vz corresponds to increasing radial dis- 
tance.) 

For our model, we tile the stellar surface in lat/long. 
To first order in SR, each tile's projected area AAp is 



AAp = {l + 2SR)n- An 



(13) 



where AQ, is the solid angle of each grid point. 

For a given orbital phase (j), we calculate SK, ST, 
and T at each grid point on the stellar hemisphere 
visible to the observer (Z > 0), along with Ro ■ ^ro. 
To include the Doppler flux variation, we assume each 
oint on the star is a blackbody - ILoeb fc Gaudil 
20031) showed that departure from blackbody emission 
changes the photometric signature of the ellipsoidal 
varia tion by only a few p ercent. We use Equation 2 
from iLoeb fc Gaudil ()2003[ ) to calculate the monochro- 
matic flux throughout the Kepler observational bandpass 
(jhttp : / /keplergo . arc . nasa. gov/ C alibr ationResponse . shtml ) 
and convolve the flux with the Kepler response function. 
(Note that the Kepler response function is given in 
wavelength space, so we had to convert it to frequenc y 
space to use the results from iLoeb fc Gaudil 120031 ) 
Using the emission calculated for each point on the 
stellar surface, we multiply each point's flux by the 
appropriate limb-darkening profile value (/(/i)//(l)) 
and sum the contributions from all grid points. Finally, 
we move to the next point in the orbit and do the 
calculation over again. Figure [2] illustrates schematically 
the appearance of the distorted star and the resulting 
photometric variations. 

Because our model is linearized in small quantities, it 
is computationally more efficient for planet-star systems 
than more general n iodels, particularly those designed 
for binary stars (e.g. lOrosz fc Hauschildtl I2000D . Where 
more general models require a few hundred thousand 



grid eleme nts to accurately model the ellipsoidal vari- 
ation (e.g. iWelsh et al.ll2010l ). our model requires only a 
few hundred for convergence to better than 0.01 ppm. To 
check the accuracy of our approximations, we compared 
the values for all calculated physical quantities (radius, 
temperature, etc.) as determined by the linearized equa- 
tions and the exact equations. For the HAT-P-7 system, 
all quantities converged to better than 10 parts per bil- 
lion. 

The ellipsoidal variation signal is convolved with the 
planet's reflection and emission in the Kepler data. 
Therefore, to flt the Kepler data, we also require a model 
for light emerging from a planet. In visible wavelengths, 
the light emerging from a close-in planet is likely dom- 
inated by reflection of stellar radiation, which suggests 
the planetary phase function should be nearly symmet- 
ric a bout (j) = .5. Although it can be more complicated 
(e.g. iCowan fc"Agol 2008) . we take a simple sinusoidal 
phase curve for the planet: 



Fp = Fo- Fi cos(2'Kq 



(14) 



where Fq is a constant term, Fi is the amplitude of the 
planet's reflected and emitted light, and both are ratioed 
to the stellar emission at mid-eclipse (when the planet is 
occulted). The sum Fq -|- Fi can be directly estimated 
from the depth of the secondary eclipse. 

Close-in planets also suffer significant tidal distortion. 
This distortion may increase the planet's phase curve at 
quadrature as a planet's projected surface area is largest 
there, a nd this effect may b e observable, particularly in 
the IR ()Cowan et al.l [201^ . However, we assume this 
effect is negligible in Kepler's bandpass. Future work 
should re- visit this assumption. 

The ellipsoidal variation depends on several key system 
parameters, although there is degeneracy between some 
parameters (/3 and q, for example). Given sufficiently 
high quality data, the ellipsoidal variation can be used 
to determine at least seven parameters: q, a, w*, sin i, 7.;, 
Kz, and (3. A fit to the planet's phase curve determines 
Fq and Fi. We expect the planet's phase curve to oscil- 
late with the orbital period and the ellipsoidal variation 
to oscillate with half the orbital period (twice an orbit). 
While they have a period equal to the orbital period, 
the Doppler fiux variations are 90° out of phase with the 
planet's phase function. Thus, in principle, analysis of 
Kepler observations should be able to distinguish these 
different components (given that the planet's phase curve 
is symmetric about </> = 0.5). 

There is no general expression relating the physical pa- 
rameters to the amplitude of the ellipsoidal variation, 
but approximating the star's shape as an ellipsoid, we 
can explicitly express the ellipsoidal variation's depen- 
dence on system parameters. Then the photometric os- 
cillations can be ex panded as a Four i er seri es. C ombin- 
ing eq uations from iMazeh fc Faigleil (|2010f) and iMorrisI 
()1985t ) gives the following series for the combined ellip- 
soidal variation, Doppler flux variations, and planet's 
phase curve ■^: 



AF 

— ;- = - Aeiiip cos(2 • 2n<j)) + Abeam sin(27r(/)) 
r 



- Arefl cos(27r0) 
where Aeiup = aeinpiqsm^ i)a^^, Abeam 



(15) 



^beam^ V c / ' 




Fig. 2. — Cartoon illustrating the phases of ellipsoidal variation. The y-axis is in arbitrary units, and the x-axis is orbital phase. The 
planet-star mass ratio is exaggerated for illustrative purposes. Phases 1 (0 = —0.25) and 3 (</> = 0.25) show the planet at quadrature, 
phase 2 ((/) = 0) shows the planet during transit (the photometric signature of which is NOT included), and phase 4 (0 = 0.75) shows the 
planetary eclipse. 



and A. 



refl — Pgeo 



-^ ) . Here pge 



is the planet's ge- 



ometric albedo, and Rp the planet's radius, aeiup de- 
pends on the gravity-darkening and limb-darkening coef- 
ficients, and abeam corrects the amplitude of the Doppler 
flux variations for shifting of flux into and out of the 
observational bandpass. Both as are of order unity 
([Mazeh & Faiglcr 2010), but more accurate estimates are 
required to provide estimates of, for example, q. More- 
over, as discussed in Section [^TTl the above Fourier series 
less accurately approximates the photometric oscillation 
for very close-in planets (a — >■ 1), as higher-order har- 
monics contribute non-negligibly. 

Although, in principle, ellipsoidal variations can con- 
strain a, i, and ji, we expect that transit observations 
(if a planet DOES transit) will provide tighter con- 
straints. Also, ellipsoidal variations are relatively insen- 
sitive to /3, and so modeling based on spectral charac- 
terization of a star may provide better estimates (e.g. 
iClaret fc^ locmcn 2011). On the other hand, when ellip- 
soidal variations can constrain q and transit observations 
sinz, the Doppler variation signal or radial velocity ob- 
servations may provide independent constraints on M^. 
Consequent to these considerations, in our analysis of 
the HAT-P-7 observations (Section 14]) , we do not fit for 
a, uJi,, i, "fi or /3 and fix these at values provided by other 
studies. 

We performed several tests to verify that our model 
works correctly. For example, in the next section, we 
compar e our model to the more genera l Wilson-Devinney 
model (|Van Hamme k Wilsonl |2007() a nd find good 
agreem ent. We also used the results from lShporer et all 
(j2011| ) to test our model for Doppler flux variations (see 



Equat ion [T2l and preceding discussion). iShporer et al.l 
(|2011l ) analyzed Doppler flux variations observed for the 
KOI-13 system and determined their amplitude to be 
9.32 ppm, corresponding to Kz = 954 m/s (see their 
Equation 1). Our model indicates that Kz = 973 m/s 
is required to produce that amplitude for t hat s ystem, 
within 2% of the result from lShporer et al.l ()2011l ). 

2.3. Comparison to Other Models 

In this section, we compare our model to previ- 
ously developed mo dels. We consid er the sinusoidal 
model proposed by iMazeh &: Faigleii (|201C0 and de- 
scribed by Equation 1151 Given the assumptions 
under which they're derived, we expect the sinu- 
soidal model to be less accurate for a ^- 1 and the 
EVIL-MC model to be less accurate a,s q — > 1. We 
also consider the publicly available Wilson-Devinney 
(W-D) model (ftp://ftp.astro.ufl.edu/pub/wilson/), 
which has a storied history and has been devel- 
oped fo r a wide variety of astro physical circum- 
stances (jVan Hamme fc Wilsonl [20071 ). Comparison to 
other models would be h elpful, but the ELC model 
(jOrosz fc Ha uschildt'l2 000|) isn't publ icly available. The 
JKTEBOP model ( Sou thworthefall [200l IS available 
(http://www.astro.keele.ac.uk/jkt/codes/jktebop.html) 
but approximates tidally distorted bodies as ellipsoids 
and so probably would not provide a more accurate 
description of tidal distortion than the sinusoidal model 
doefl 

^ After this paper was accepted for publication, Dr. Jan Buda j 
made us aware of another relevant model described in'BudaJl 1)20111 ) 
and available at http://www.ta3.sk/~budaj/shellspec.html 



TABLE 1 

Model fit parameters from our analysis 



par am. 


value (fixed D) 


value (var. D) 


9 


(1.10 ±0.06) X 10-^ 


(0.99 ±0.07) X lO--^* 


D 


61 ± 3 ppm 


65 ± 2 ppm 


-^ day 


2680+^^ K 


2700 ± 10 K 


Fi 


30 ± 1 ppm 


32 ± 1 ppm 


t'o - Fi 


or 1 ppm 


3 ± 3 ppm 


Kz 


300 ± 70 m/s 


300 ± 70 m/s 



Note. — The middle column shows best-fit values for fixed D 
and fixed/ variable Kz, while the rightmost column shows values 
for variable D. 

Although the W-D model is widely applicable, for 
modeling planet-induced ellipsoidal variations, its nu- 
merical precision is limited to a few tens of ppm (R. E. 
Wilson, private communication, 2012). Consequently, in 
the comparison below, the smallest g-value we consider 
is 0.05, which corresponds, for example, to a 50 Jupiter 
mass body orbiting a solar-mass star. For the range of 
relevant a-values, g-values more appropriate to planets 
{q < 10"'^) produce ellipsoidal variations below the W-D 
model's numerical precision. In any case, the range of q 
available is sufficient for our purposes. 

For the comparison, we fix several model parame- 
ters. Neither the W-D nor the sinusoidal model al- 
low quadratic limb-darkening, so we assume linear limb- 
darkenin g for the comparison, with a coefficient u = 
0.551 ('Cl aret fc Blo cmcn 2011). (For EVIL-MC, this as- 
sumption is equivalent to 71 = 0.551,72 — 0.) We take 
the gravity-darkening coefiicient to he /3 — 0.071 (cor- 
responding to 5 = 4/3 = 0.284 for W-D). With these 
parameters, aeiup — 0.15 (15 + w) (1 + g) / (3 — w) = 
1.223. We do not include Doppler flux variations and 
reflected/emitted light from the planet for this compari- 
son. For the sinusoidal model, this assumption requires 
Abeam = Arefi = 0. For the W-D modcl, we set the 
secondary's luminosity (L2) to zero. We also assume no 
stellar rotation. Unless specified below, all other system 
parameters are fixed at the values in Table [2l 

First, we compare results from the three models for 
a range of q and a = 3, as illustrated in Figure [3] (a). 
To determine the overall normalization for the sinusoidal 
model, we added an offset value to Eg uationlTCland use d 
a Levenberg-Marquadt (LM) scheme ()Markwardtll2009[ PI 
to find the value that provided the best agreement be- 
tween the sinusoidal and W-D models. 

As illustrated in Figure [3] (a) , agreement between the 
W-D and EVIL-MC models is better than 2.1% of the 
total ellipsoidal variation for all q illustrated, even though 
the EVIL-MC model is derived under the assumption of 
small q. For q — 0.05, the two models agree to 1.1% of 
the ellipsoidal variation, corresponding to a difference of 
about 40 ppm. This discrepancy is near the numerical 
precision limit of the W-D model and so is as good as the 
agreement can be. These results indicate the EVIL-MC 
model is sufficiently accurate to model tidal distortions 
even in binary systems with stars of comparable mass. 
By contrast, the sinusoidal model agrees with the W-D 
model to only about 10% of the total variation. 

Next, we compare results for a range of a and q = 0.05, 
as illustrated in Figure |3] (b) . Agreement between the 

^ We used Craig Markwardt's mpfit.pro IDL routine, available 
at |http://www. physics. wisc.edu/~craigm/idl/fitting. html 
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Fig. 3. — Ellipsoidal variations predicted by the Wilson-Dcvinney 
(W-D) (black lines), the EVIL-MC (blue), and the sinusoidal 
(Equation I15|l (red) models for a range of q and a. The models 
here do NOT include Doppler fiux variations, refiected/emitted 
light from the secondary (planet), or the transits/eclipses. (The 
transit phase occurs to the outside of the black, vertical lines.) 
The models DO include limb- and gravity-darkening, (a) Predic- 
tions for fixed a = 3 and q equal to 0.05 (dash-dot-dot-dot lines), 
0.1, (dash-dot), 0.5 (dash), and 1 (solid), (b) Predictions for fixed 
q = 0.05 and a equal to 2 (solid lines), 3 (dash), and 5 (dash-dot). 
The difference between the EVIL-MC and W-D models is 2.1% or 
less of the total variation and is as good as the agreement can be, 
given limits on the W-D model's numerical precision. By contrast, 
the difference between the sinusoidal and W-D models is usually 
greater than 10% of the total variation. 

EVIL-MC and W-D models is better than 1.5% of the 
total variation, while agreement between the W-D and 
sinusoidal model is no better than 8% and as bad as 20% 
(for a = 2). As expected, the sinusoidal model is less 
accurate as a — >■ 1 as higher order Fourier components 
contribute more. Whether there exist planets with a — 2 
and the sinusoidal model can be applied to them remains 
to be s een (tidal decay of th eir orbits wou ld pro bably be 
rapid - iLevrard et al.li2009l : IJackson et al.ll2009D . but the 
Kepler mission has announced candidates with a ~ 2. 

We can ask how accurate are estimates of system 
parameters from the sinusoidal model, particularly the 
mass ratio. Figure |4] illustrates the accuracy of the q- 
value estimated using the sinusoidal model. For that 
figure, we calculated ellipsoidal variations for a range of 
a- and g-values using the EVIL-MC model (again, ne- 
glecting Doppler fiux variations or reflected/emitted light 
from the planet). Then, we used an LM scheme to deter- 
mine a best-fit Aeiup (and offset value) for each modeled 
ellipsoidal variation and estimated q from Aeiup by using 
the assumed values for all other system parameters (a, 
aeiup, etc.). 

Figure|4]shows the ratio of the g- value estimated in this 
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Fig. 4. — Rati o o f the g-value estimated using the sinusoidal 
model (Equation HSU and the actual g- value (shown along the x- 
axis) for a range of q and a. The sinusoidal model always underes- 
timates the actual mass ratio by a few percent, and the estimate's 
accuracy degrades for a — > 1 as higher order Fourier components 
contribute more. 

way to the actual value. The sinusoidal model typically 
underestimates g by a few percent, and, as expected, the 
estimates become less accurate for small a. Although es- 
timates of q from, for example, Kepler data are likely to 
be less accurate than a few percent, estimates of q using 
the sinusoidal model may be systematically smaller than 
the actual q-values. Depending on how the modeling is 
done, inaccuracies in the estimation of q may cause es- 
timates of other system parameters to be systematically 
inaccurate as well. In any case, we confirm that the sinu- 
soidal model should generally be sufficiently accurate to 
distinguish planetary companions from low mass stellar 
companions. 

3. OBSERVATIONS OF THE HAT-P-7 SYSTEM 

The HAT-P-7 planetary system was discovered by the 
HATNet survey and was the second planet discovered 
in the Kepler field of view. The system is composed of 
an F-type star (IVh = 1A7Mq, R^, = IMRq) and a 
gas giant planet {Mp = 1.78M./„,p, R^ = 1.36Rjup) in 
a 2.2 day circular orbit ()Pal et al.l [20080 . In the Ke- 
pler bandpass, the star has a magnitud e K^, = 10.5, 
relativ ely bright among Kepler targets. iBorucki et al.l 
(j2009| ) analyzed the first ten days of Kepler data (quarter 
0, QO), detected a secondary eclipse depth of 130 ±11 
ppm, and estimated a d ayside temperature of 2650 K. 
iChristiansen et all ()2010[ ) analyzed observations of HAT- 
P-7 b's secondary eclipse from the EPOXI mission and 
put upper limits on its depth at 0.055%. They also an- 
alyzed Spitzer secondary eclipses taken throughout the 
IR and found brightness temperatures in the different 
b andpasses from 225 K to 3190 K. 

iWelsh et al.l ()2010D discovered ellipsoidal variations in 
the Kepler Ql data, with an amplitude of 37.3 ppm. 
They also estimated the planet's phase curve has an am- 
plitude of 31.9 ppm and day- and nightside tempera- 
tures of 2885 and 2570 K, resp e ctively . The discrep- 
ancy between the IBorucki et al.l (|2009l ) result and the 
IWelsh et al.l ()2010[ ) result may arise from the considera- 
tion of ellipsoidal variation in the latter analysis and/or 
the inclusion of more data (Ql data span 30 days, as 
compared to QO's 10 days). 
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Fig. 5. — (a) Kepler observations of HAT-P-7 system from Q2 
with outliers filtered out (see text), (b) Observations from QO-2 
phased together and binned to 30- minute bins (Xs). Our best fit 
model curve (solid line) is also shown, the best-fit ellipsoidal vari- 
ation is shown as a dashed curve, the planetary phase curve as 
a dash-dot line, and the Doppler fiux variations (with Kz = 300 
m/s) as the dash-dot-dot-dot line. The planet's eclipse is high- 
lighted in grey and is not fit by our model. Best fit parameters 
are shown in Table [l] Our originally estimated uncertainties (~4 
ppm) are re-scaled from the initial values by 1.8, the square root 
of the best-fit reduced x = 3.25, to give uncertainties ~8 ppm. 
(c) Residuals between the best-fit model and the data are nearly 
normally distributed about 0. 



I Winn et al.l ()2009al ) conducted Rossiter-McLaughlin 
observations of HAT-P-7 and found the planet is in a 
near polar or even retrograde orbit about its star, with 
an angle between the stellar rotation and orbit normal 
vectors projected onto the sky- plane of 182.5±9.4°. They 
also pointed out that the unusually low projected rota- 
tional velocity of the star for its type and age suggests we 
are observing the star nearly pole-on. (Note that even for 
the e stimated deproject ed rotation velocity of ~15 km/s 
from lWinn et al.ll2009al HAT-P-7 is stih a "slow-rotator" 
for the purposes of our linearized model.) 

For our analysis, we obtained the short-cadence (1- 
minute observing cadence) Pre-Search Data Condi- 
tioned (PDC) light c urves from the MAST archive 
(http://archive.stsci.edu/kepler/), and we analyzed the 
QO-2 data. (I'he other data publicly available at the time 
of our analysis from Q3 exhibited more complex system- 
atic trends, so we did not include them in our analysis.) 
We first removed outlying photometric points that were 



explicitly flagged in the data files as anomalous by the 
Kepler team. We then binned the remaining data in 30- 
minute bins, calculated each bin's standard deviation, 
and threw out data points more than 4 standard devia- 
tions from the mean in each bin. The Q2 data filtered 
in these ways are shown in Figure [5] (a). These data 
still clearly exhibit both long-term trends and correlated 
noise. 

We attempted to remove these trends. First, for the 
QO data, we masked out all the transits (5 transits) and 
fit a 4th-order polynomial to the remaining data. (3rd 
and 5th-order polynomials gave equivalent results within 
uncertainties.) We calculated the standard deviation a of 
residuals between data and the trend curve and dropped 
points that lay more than 4-(t from the trend curve. 
We re-fit a 4th-order polynomial to these screened data 
and iterated this procedure until all A-a outliers were re- 
moved. Then, we divided the data (including transits) 
by the final trend curve. We performed the same de- 
trending for the Ql and Q2 data. Analyzing these three 
quarters together nearly qua druples the number of orbits 
examined over the analysis of lWelsh et al.l (|201Clf ) and sig- 
nificantly improves the accuracy of the estimated system 
parameters. 

After detrending the data, we phased and stacked 
them, assuining a n orbital period of 2.204733 days 
(jWelshet al.ll2010[ ). We then binned the data into 30- 
minute wide bins, determined a median for each bin, and 
took 1.4826 X the median absolute de viation (MAD) a s 
the standard deviation for each bin (jBevingtonl Il969[ ). 
We then threw out points in each bin more than 4-cr 
from the median. We then took the median of the re- 
maining data in each bin. For the uncertainties, we took 
1.4826 X MAD divided by the square root of the number 
of points in each bin - such uncertainties were typically 
4 ppm. However, sys tematic trends or correlated noise 
still pervade the data ()Pont et al.ll2006D , producing scat- 
ter larger than 4 ppm. 

To estimate the size of this scatter, we determined an 
initial best-fit model usin g a Levenberg-Marquadt algo- 
rithm (jMarkwardtl 120091 1. which gave a reduced x^ = 
3.25, indicating the scatter was indeed underestimated. 
We re-scaled the error bars by \/3.25 = 1.8, giving un- 
certainties ^8 ppm. 

Finally, we calculated the overall normalization of the 
data by taking the mean of the data during the eclipse 
phase, when only the star is contributing flux. (Es- 
timated variation of the system brightness during this 
phase is less than 0.1 ppm, and so variations in these 
data are dominated by intrinsic scatter.) We divided all 
the data through by this value (for Figure [5] (b), we sub- 
tracted 1.0 from the data). These are the final data we 
analyzed and are shown in Figure [5] (b) (with the tran- 
sit near phase and 1 masked out), along with our best 
model curve (see below). The contributions from ellip- 
soidal variations, the planetary phase curve, and Doppler 
flux variations are also shown. 

4. ANALYSIS 

We conducted a suite of Markov-chain Monte Ca rlo 
(MCMC) analyses, using Gibbs sampling (|Fordll2005D to 
fit the model parameters, q, Kz^ -Fq, and Fi (Table [T]). 
In some of the model fitting, though, we also held Kz 
constant, and the sum Fq + Fi was constrained by the 
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12.85" 
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'1\ 


6350 K" 


m 


0.26'" 


log(g) 
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I 
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p 


2.204733 days'' 


(71.72) 


(0.314709, 0.312125)*^ 
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iWelshcFal? ('2013) 



tPal ct al. ( 2008) 
'^Dete rmined from 
IClaret 



interpolation 
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among the values 



eclipse depth (see below). We held all other parameters 
fixed for all modeling (Table [5]). 

The eclipse depth provides a constraint on the max- 
imum of the planetary phase curve D — Fq + Fi. We 
estimated the eclipse depth by fitting a straight line be- 
tween the points on either side of the eclipse. Then, 
we took the eclipse depth to be the difference between 
that value and 1.0 (the normalized value during eclipse), 
giving D = 61 ± 3 ppm. (The difference between the 
actual maximum in the planetary phase curve and the 
maximum estimated this way is considerably less than 
the scatter in the data.) We conducted two sequences 
of MCMC analyses: (1) with D held constant (best-fit 
parameters for which are in the first column in Table [TJ 
and (2) allowing D to float but with a x^-penalty for 
departures from 61 ppm (second column). (For the lat- 
ter analysis, Kz was allowed to float as well.) Compari- 
son of the two sequences below highlights the degeneracy 
between constraints on the ellipsoidal variation and the 
planetary phase curve, but we focus our discussion on 
the analyses with D — const since they give a g-value 
consistent with previous analyses. 

In theory, the ellipsoidal variation signal depends on 
the relative orientation of a;^, and the orbit normal vec- 
tor. However, we tried different relative orientations and 
found the data cannot distinguish between an a;^ that 
points directly at the observer (parallel to Z - Figure [1]) 
and any other orie ntation allowed by other constraints 
(jWinn et al.l l2009al ) . Thus, we assumed the w^ vector 
points directly at the observer in our modeling (w^^HZ). 

For the MCMC fitting, we used 5 Markov chains, each 
with 5 X 10"^ links, which we show below suffices for good 
convergence of the model parameters. (We also con- 
ducted MCMC analyses with 5 x 10"* links which con- 
firmed the shorter chains had converged.) For each jump 
transition, we chose at random either 1 or 2 parameters 
to vary. We discarded the first 20% of each chain. Oth- 
erwise, the resulting distributions of best-fit parameters 
might be skewed by our initial choice of parameter values. 
For sampling the paramet er space, we took the Gaussian 
distribution suggested by iFordl (|2005t ) for the candidate 
transition probability, with a width /3 for each parameter 
such that the fractio n of accept ed transitions was ^0.25. 
(See Equation 12 in Hb^ [20051 ) 

We also checked that our analysis technique can ac- 
curately recover system parameters by generating sev- 
eral synthetic data sets designed to mimic the raw Ke- 
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Fig. 6. — The mean values for each parameter from each of the 
5 Markov Chains as a function of link number (after the first 20% 
of each chain was dropped). Each linestyle represents a different 
chain. The left column ((a)-(c)) showrs the chains for wrhich Kz 
is taken as a free parameter, while the right column ((d) & (e)) 
shows the chains for which Kz = 213.5 m/s. The vertical lines in 
each panel aligned with x = 3,500 represent the smallest standard 
deviation from among all the chains, and the means for all chains 
converge to within those deviations. 

pier data with ellipsoidal variations, planetary emis- 
sion/reflection, Doppler flux variations, and the same 
gaps in time and scatter. We consistently recovered the 
assumed system parameters when they were recoverable 
(see discussion of Kz estimate below) . 

For the sequence of analyses with the eclipse depth D 
held constant. Figure [5] illustrates the convergence of the 
mean for each of the 5 chains. For each parameter, the 
distribution from each chain provides a slightly differ- 
ent mean value, but the differences between the various 
mean values are all smaller than the smallest standard 
deviation for any one chain. For example, in Figure|n](a), 
the largest difference in the final mean q-value between 
different chains is 5.3 x 10~^, while the smallest stan- 
dard deviation from among all the chains is more than 
10 times larger, indicating the chains have all converged 
within uncertainties. 

Figure [7] shows the distributions of best-fit parame- 
ters for the two sequences of model-fitting with con- 
stant Z?, one with Kz variab l e, th e other with Kz 
fixed at 213.5 m/s ()Pal et al.l I2008D . The mean of 
each distribution is taken as the best-fit value, and the 
standard deviation is the uncertainty. Our best-fit q 
((1.10 ± 0.06) X 10 ~^) is smalle r but c onsistent (within 
2-a) with that of IWelsh et all (poTol ) (1.190 x lO^^). 
Assuming M^ = 1A7Mq, our g- value corresponds to 
M„ = 1.62Mj„„. 

iPal et al.l (|2008[ ) estimated Kz = 213.5 m/s, but our 
best-fit value is 300 ± 70 m/s. This latter value corre- 
sponds to Doppler flux variations of only about 4 ppm. 
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Fig. 7. — The distributions of the model fit parameters result- 
ing from our MCMC analysis, with Kz variable (solid) and fixe d 
(dashed), (a) Mass ratio, q - For A/* = 1A7Mq dPal et al.1 12008'). 
our best-fit q gives Mp = 1.62Af j^^. (b) Amplitude of the plan- 
etary phase curve, Fi - The sum Fq + Fi, the emission from the 
planet's dayside hemisphere, is held fixed at the estimated eclipse 
depth, 61 ppm. With the best-fit Fi illustrated, _Fo — Fi, the 
emission from the planet's nightside hemisphere, is nearly 0. (c) 
Amplitude of the stellar reflex velocity, Kz - Even for Kz = 300 
m/s, the signal from the Doppler flux variations has an amplitude 
of only about 4 ppm, twice as small as the intrinsic scatter in the 
data. Thus, the best-flt Kz is very sensitive to the scatter and has 
large uncertainties. 

below the intrinsic scatter in the data. We conducted 
numerical tests to see whether we could, indeed, have re- 
covered a Kz = 213.5 m/s and found that we could only 
recover Kz for scatter less than or comparable to the 
Doppler signal. IWelsh et al.) ()2010t ) estimated Kz = 212 
m/s by steering Kz toward 213.5 m/s via a x^ penalty 
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for deviations (although in Figure 3 of that study, the 
asymmetry that should result from the Doppler signal 
seems absent). In any case, as illustrated in Figure [7] 
(a) and (b), the best- fit values for other parameters are 
insensitive to Kz- 

As discussed above, our estimated eclipse depth for the 
planet is 61 ± 3 p pm, as compared t o the eclipse depth 
of 85.8 ppm from iWelsh et al.l ()2010[) . This discrepancy 
arises from our use of more data than used in that pre- 
vious study. A preliminary analysis of Q l data alone 
yielde d an eclipse depth similar to that of I Welsh et al.l 
(J201C1( ). Assuming the planet's day side emits as a uni- 
form blackbody, our depth corresponds to a dayside tem- 
perature of 2680^20 K' which is almost 20 K smaller 
than the average dayside temperature from IWel sh et all 
y^QlO). The disagreement with the eclipse depth from 
iBorucki et al.l (|2009l ) probably arises for similar r easons. 
As a further confirmation of our estimate, iMislis et al.l 
(|2012D conducted an analysis of some of the same data 
as we and found a similar eclipse depth (see their Figure 

For our analyses with fixed D, we found that the night- 
side emission Jn — fi — 0, as com pared to the 22.1 ppm 
estimated bv iWelsh et al.l ()2010[ ). Partly, this disagree- 
ment is due to the fact tha t we do not explic itly ana- 
lyze the transit phase, while IWelsh et al.l (J2O1O0 do, and 
partly, it is due to our choice of planetary phase function: 
IWelsh et al.l (|2010f ) chose a planetary emission/refiection 
relationship that produces a shallower drop off in plan- 
etary flux than our function as departs from 0.5. 
Both estimates for the nightside emission are model- 
dependent, though. 

FigureJHJillustrates the results of the MCMC analysis in 
which D was allowed to float. (Note: convergence of the 
model parameters for this analysis required 5 x 10** links 
in each chain.) Our MCMC analysis drives q to smaller 
values and D to larger values than when D is held flxed. 
Unfortunately, because we re-scaled our uncertainties to 
force x^ ~ Ij we cannot use the x^-values from the dif- 
ferent MCMC sequences to determine whether letting 
D float provides a statistically better model flt. How- 
ever, the fact that the g- value for the sequence with 
fixed D more closely matches previous constraints sug- 
gests that is the more appropriate model. In any case. 
Table [I] shows the planetary parameters corresponding 
to the best-fit values for variable D. 

We also conducted numerical experiments for which we 
created synthetic datasets with the same best-fit param- 
eters produced by the previous MCMC analysis (column 
1 in Table ll]). We added Gaussian noise to these syn- 
thetic datasets, with a scatter of 8 ppm. We applied the 
same MCMC analysis in which we allowed D to fioat, 
and the analysis would often drive q to smaller values 
and D to larger values than assumed, depending on ex- 
actly where the noisy data points ended up. These re- 
sults highlight the degeneracy between the best-fit q and 
planetary phase function and show that it can depend 
sensitively on the scatter in the data. 

The dependence of the derived g-value on the as- 
sumed plane t ary p hase function has been considered by 
IMislis et al.l (|2012| ). For any planetary phase function 
that is symmetric about (p = 0-5, there will necessarily 
be some degeneracy between the solution for the phase 
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Fig. 8. — The distributions of the model fit parameters resulting 
from our MCMC analysis, with Kz variable. The eclipse depth 
D is also allowed to float but with a x^-penalty for departures 
from 61 ppm. (a) Mass ratio, q - The MCMC analysis drives q to 
smaller values when D is allowed to float, producing a best-fit q = 
0.99 ± 0.07 X 10^-^. (b) Amplitude of the planetary phase curve, Fi 
- This value also goes up as D increases, producing a best-fit Fi = 
32 ± 1 ppm. (c) Amplitude of the stellar reflex velocity, Kz - This 
parameter is essentially unchanged and has a best-fit value Kz = 
300 ± 70 m/s. (d) Eclipse depth D - The MCMC routine drives 
this parameter to larger values than our best estimate, producing 
a best-fit value Z) = 65 it 2 ppm. 
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Fig. 9. — The values for mass ratio q and the amplitude of the 
planetary phase function Fi sampled during the MCMC analysis. 
The strong positive correlation between these parameters arises be- 
cause an increase in flux from a planet-star system when the planet 
is near quadrature can be attributed to increasing i<b /reducing Fi 
(their sum is constrained) and reducing q or vice versa. 

function and q. For example, a model fit to light emerg- 
ing from a planet-star system exhibiting ellipsoidal vari- 
ations can enhance the peaks near (jj = 0.25 and 0.75 
by increasing the baseline planetary flux (Fq), subject to 
constraints on the eclipse depth, or by increasing q. Ad- 
ditional constraints on the planetary emission from other 
phases cannot completely remove this degeneracy. Addi- 
tional degeneracies between, for example, the planetary 
emission and the transit parameters should emerge from 
analysis of the transit phase. 

Even for D ~ const., this degeneracy persists. Figure 
[3] illustrates the degeneracy between q and Fi and shows 
the values sampled by the MCMC routine (with variable 
Kz) when D is held constant. Since D = Fq + Fi — 
const., an increase in the signal from the system when 
the planet is near quadrature can be attributed to either 
increased q or increased Fq /decreased Fi. A similar de- 
generacy does not exist for Kz since the Doppler flux 
variations aren't symmetric about </> = 0.5. Taken alto- 
gether, our results show that, while Kepler and CoRoT 
data may provide important constraints on planetary 
albedo and energy budget, constraining these properties 
requires including ellipsoidal variations. The different 
contributions cannot be completely disentangled. Fu- 
ture work should consider more completely the influence 
of alternative planetary phase functions. 

5. DISCUSSION AND CONCLUSIONS 

We have developed a new model for ellipsoidal varia- 
tions induced by close-in planets, the EVIL-MC model. 
EVIL-MC employs several approximations suited for 
planet-star systems and is thus computationally more 
efficient than other more general models and more ac- 
curate than simpler, semi-analytic models. For example, 
the W-D model takes about 0.5 s to run each of the exam- 
ples in Section [231 while our EVIL-MC model (in IDL) 
runs in less than 0.05 s. Also, after our HAT-P-7 data 
were detrended and binned, we performed the entire suite 
of MCMC calculations (25,000 evaluations of the EVIL- 
MC model) described in Section H] in about 20 minutes. 
This increase in efficiency makes our model well-suited 
for analyzing the mountain of Kepler and CoRoT data 



still pouring in. 

The EVIL-MC model has some important limitations. 
It is not designed for systems with a mass ratio q ^ 1 
since tidal distortions are not small for those systems, al- 
though our comparison to the more general W-D model 
shows agreement at about 2% even for large q-values 
(Section [131). 

The EVIL-MC model may not be sufficiently accurate 
for rapidly rotating stars where the rotational oblateness 
is large. Rotationally induced gravity darkening at the 
equators of such stars may imprint a discernible signature 
on the transit light curves of companion planets, analysis 
of which can reveal the misalignment between a planet's 
orbit and the stellar equator. Such analyses have been 
conducted for the KOI-13.01 system (jBarnes et al.ll201ll : 
ISzabo et al.l[20Tlll . 

EVIL-MC also does not currently include the transit 
and eclipse phases for a planetary system, and a future 
version will also include these phases. However, a prelim- 
inary analysis shows that tidal distortion has a negligible 
(< 0.1 ppm) influence on the transit light curves for typ- 
ical planetary systems. 

Accurate determination of planetary phase curves from 
Kepler and CoRoT data requires consideration of the 
ellipsoidal variations, and there can be degeneracies be- 
tween the contributions from the ellipsoidal variation and 
the planetary phase curve. Planetary phase curves are 
diagnostic of atmospheric temperatures and dynamics, 
and a complex story of coupled chemistry, dynamics, and 
radiation is emerging, motivated largely by IR observa- 
tions of_2la£etarY_P,li2;Se curves and eclipse depths (see, 
e.g.. iKnutson et aLll2Q10D . Results from the Kepler and 
CoRoT missions will add to this picture and, when com- 
bined with Spitzer observations, will give a much fuller 
picture of the atmospheric energy budgets of close-in 
planets. 

From our analysis, we can draw some tentative con- 
clusions regarding HAT-P-7 b's atmosphere. Given its 
proximity to its host star, the planet is probably tidally 
locke d, and the same side of the planet always faces the 
star (j Jackson et al.l [20081 ). Consequently, atmospheric 
circulation is required to transport stellar heating from 
the day- to the nightside. Our estimated dayside emis- 
sion 61 ppm corresponds to a brightness temperature of 
2680 K. Our estimated minimum for the planet's night- 
side emission Fq — Fi lies below the sensitivity of our 
analysis, ~ 4 ppm, suggesting the nightside brightness 
temperature in the Kepler band is less than 1970 K. 

This result might indicate much of the stellar heat- 
ing on the dayside is radiated to space before it can be 
transported to the nightside. This result is also qual- 
itatively c onsistent with mode ls of the hottest close- 
in planets (jFortnev et al.l [20081 ) and with analyses that 
sugge st HAT-P-7 b has an atm ospheric thermal inver- 
sion ([Christiansen et all 120100 . which is often corre- 
lated with a high atmospheric temperature for close- 
in planets ([Knutson et al.l [2010[ ). However, determin- 
ing the precise implications of this result for the at- 
mospheric circulation requires detailed modeling. It is 
worth noting that the day-night brightness temperature 
contrast inferred here ( 710 K) is similar to that in- 
ferred for WASP-12 b (iCowan et al.l I2012D but greater 
than those of cooler hot Jupiters, inc luding HP 189 733 b, 
^Knutson et al.,,20071 HD 209458b ([Cowan et al.|[2007.) . 
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and HD 149026b (jKnutson et all I2OO90 (although the 
latter has an error bar that does allow relatively large 
values). These measurements are not all at the same 
wavelength, which complicates the interpretation. An 
additional complication is that our eclipse depth may 
also include contributions from atmospheric scattering 
of light by clouds, although estimat ed optical albedos 
of hot Jupiters are highly uncertain (jRowe et al.l I2Q08I : 
ICowanfc Agol[20Tll) . 

To help place our results regarding HAT-P-7 b in 
context, we ran some preliminary dyna mical+radiative 
calcu lations using the SPARC model ([Showman et al.l 
l2009f) . The HAT-P-7b model atmospheres were con- 
structed assuming a solar metallicity atmosphere in ther- 
mochemical equilibrium for cases with and without TiO, 
which can absorb stellar radiation hig h in the atmosphere 
and produce a temperature inversion (jChristiansen et al.l 
[20im . 

For the model with TiO in the atmosphere, stellar 
heating is deposited higher in the atmosphere, where the 
timescale for radiation of stellar heating to space is rela- 
tively short, and consequently the model predicts a large 
day-night contrast: a day side emission of 74 ppm and a 
nightside emission of only 1 ppm. For the model with- 
out TiO, stellar heating is deposited deeper in the atmo- 
sphere, where the radiative timescale is relatively long, 
and so the model predicts a smaller day- night contrast: 
a dayside emission of 52 ppm and nightside emission of 9 
ppm. Comparison with our observations suggests HAT- 
P-7 b's real atmosphere might occupy a point in param- 
eter space somewhere between these models. Our results 
here suggest there are still important, unanswered ques- 
tions about HAT-P-7 b. The planet is one of the hottest 
hot Jupiters known and orbits one of the brightest Ke- 
pler targets, and so further study of the planet may prove 
particularly useful for understanding hot Jupiter atmo- 
spheres. 

Kepler and CoRoT observations will provide numer- 
ous opportunities for similar phase curve analyses. The 
closer a planet to its host star, the more stellar radia- 
tion it will receive, probably leading to greater reflection 
and/or thermal emission. The tidal distortion and el- 
lipsoidal variation of its host star could also be larger. 
Accurate determination of phase curves for the closest- 



in planets will therefore require inclusion of the stellar 
ellipsoidal variation. Phase curve fitting without it may 
produce erroneous results. By contrast, ellipsoidal varia- 
tion of a star has less influence on determination of plan- 
etary phase curves from Spitzer observations because the 
planet-star contrast for most extrasolar systems is much 
larger in the IR. 

Ellipsoidal variation analysis may provide other key 
information about e xtrasolar systems. For example, 
ILoeb fc Gaudil (j200^ first suggested Doppler flux vari- 
ations would be an important source of variability for 
Kepler observations. Equation [T^] shows that Doppler 
variations (or radial velocity observations) has a different 
dependence on the system parameters than ellipsoidal 
variations. Potentially, transit observations would give 
the orbital period P and orbital inclination sini, ellip- 
soidal variations would give the planet-star mass ratio q, 
leaving only the stellar mass unknown in Equation ll2l 

Single planets close enough to their star to induce mea- 
surable ellipsoidal variations are likely to have negligible 
orbital eccentricities. However, if primordial eccentrici- 
ties remain or interactions with other planets keep eccen- 
tricities non-zero (and the planet's semi-major axis isn't 
aligned along the line of sight), the planet-star orbital 
separation will be different at each quadrature. Conse- 
quently, the ellipsoidal variations at one quadrature may 
exceed that at the opposite quadrature, and the differ- 
ence may help constrain the orbital orientation and ec- 
centricity (iMislis et al.ll2012| ). 

Moreover, given the number of planets likely to be dis- 
covered by the Kepler and CoRoT missions, follow-up 
resources to determine the system parameters will be 
limited, so the ability to determine some of the parame- 
ters from mission photometry alone will be a tremendous 
boon. Thus, ellipsoidal variation analysis of Kepler and 
CoRoT systems promises to reveal a unique wealth of 
information. 
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